Particles on curved surfaces - a dynamic approach by a phase field crystal model 
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We present a dynamic model to study ordering of particles on arbitrary curved surfaces. Thereby 
the particles are represented as maxima in a density field and a surface partial differential equation 
for the density field is solved to the minimal energy configuration. We study annihilation of disloca- 
tions within the ordered system and premelting along grain-boundary scars. The obtained minimal 
energy configurations on a sphere are compared with existing results and scaling laws are computed 
for the number of excess dislocations as a function of system size. 

PACS numbers: 



Problems related to optimal ordering of particles on 
curved surfaces date back to the classical Thomson prob- 
lem pQ to find the ground state of N particles on a sphere 
interacting with a Coulomb potential. A classic theorem 
of Euler shows for a triangulation of the surface in which 
nearest neighbors are connected, that 5^(6 — i)vi — 6x, 
with Vi as the number of vertices with i nearest neighbors 
and x as the Euler characteristic of the surface. Thus for 
surfaces with the topology of a sphere (% = 2), besides 
the expected triangular lattice with six-fold coordination, 
which would give the optimal packing in a plane, there 
must be at least 12 five-fold disclinations present. With 
each disclination an extra energy is associated (relative 
to a perfect triangular lattice in flat space) which grows 
proportional to r 2 , with r as the radius of the sphere. 
For a fixed lattice constant a we have N ~ (r/a) 2 . Thus 
for large N mechanisms are expected which reduce this 
extra energy by changing the ground-state configuration. 
One mechanism is a buckling transition of the disclina- 
tions, which form sharp corners and turn the sphere into 
a polygon. The transition depends on Young's modu- 
lus Y and bending rigidity b via the Foppl-von Karman 
number Yr 2 jb and is intensively studied for viral cap- 
sids where protein subunits play the role of the particles 
[2j [3] . In cases where large surface tension limits signifi- 
cant buckling the energy can be reduced by introducing 
grain-boundary scars. Realizations are for example water 
droplets in oil, which are coated with colloidal particles 
[4] . Such coated droplets are potential drug delivery vehi- 
cles [5j[6]. Similar configurations occur if a jammed layer 
of colloidal particles separats two immiscible fluids form- 
ing a so-called bijel [7], which has potential applications 
as an efficient micro-reacting media. A large number of 
ordered particles on curved surfaces is also required for 
fabrication of nanostructures on pliable substrates, e.g. 
to make foldable electronic devices [8]. For all these ap- 
plications a detailed understanding of the grain bound- 
ary scars is of interest as they may be sources of leaks, 
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influence mechanical properties or lead to failure in elec- 
tronic devices. We introduce an efficient way to compute 
the dynamics of these grain boundary scars and disloca- 
tions associated with them and provide an approach to 
compute optimal ordering of many particles on arbitrar- 
ily curved surfaces. As the grain boundary scars belong 
to the thermal and mechanical equilibrium our approach 
is based on energy minimization with the geometric frus- 
tration resulting from the curved surface incorporated. 

For 2 < N < 100 there is agreement of all numeri- 
cal and theoretical methods for the Thomson problem, 
suggesting that the global minimum configuration has 
been found. However, for large N, owing to an expo- 
nential growth in local minima [9], finding global min- 
ima becomes extremely difficult. Grain boundary scars 
are expected for N > 360. Numerical approaches to 
solve such problems are typically based on genetic al- 
gorithms, steepest decent minimization or coarse grained 
approaches, in which the elasticity field between grain 
boundary scars is solved [TO] [11]. All approaches are 
devoted to finding the ground state. Dynamic models 
have been considered [12] which allow us to describe ex- 
perimentally observed dislocation glide within the grain 
boundary scars p~3j [14] . We will introduce an approach 
without any coarse graining by directly addressing the 
dynamic evolution and rearrangement of the particles on 
an arbitrary curved surface. Our approach is based on a 
free energy functional for a number density. In the plane 
such free energy functionals have been used to character- 
ize patterns. The simplest possible form of a free energy 
which produces periodic structures in a domain ft reads 

m = J^-\v P \ 2 + l\Ap\ 2 + f(p)dn (i) 

with p as the number density and f{p) = \(l — e)p 2j r\p^ 
as a potential with a parameter e. The equilibrium state 
for Q = M 2 has a perfect six-fold symmetry. Evolutional 
laws associated with this energy are the L 2 -gradient flow 
dtp = —dT/Sp, the Swift-Hohenberg model [15 , and 
the H~ x gradient flow dtp = AST/dp, the phase field 
crystal (PFC) model [16] . The evolutions naturally con- 
tain elastic energy, as an expansion of the free energy 
around the equilibrium period spacing results in the po- 
tential energy of a spring, ieHooke's law. As the energy 
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is rotationally invariant arbitrary orientationtions of pe- 
riodic structures can emerge. Furthermore the model 
allows the formation of dislocations, which occur when 
two periodic structures of different orientation collide or 
when it is energetically favorable for them to nucleate. 
The PFC model has been used to simulate various crys- 
tal growth phenomena including epitaxial growth, nu- 
cleation, commensurate-incommensurate transitions and 
plastic deformations. In [17] is shown how the model can 
be derived from a microscopic Smoluchowski equation via 
dynamical density functional theory. Formulating the en- 
ergy m Eq. on a curved surface T leads to 

? T \p] = J "IVrH 2 + \\&rp\ 2 + f(p) dT (2) 

with p as the number density on T, as well as Vr and 
Ar as the surface gradient and surface Laplacian, respec- 
tively. We will use the H -1 gradient flow of this energy 
to solve the generalized Thomson problem and to analyze 
the dynamics of rearrangements of particles on a curved 
surface. The equation, written as a system of three sec- 
ond order equations, reads as 

d t p = A r u (3) 
u = A r v + 2v + f'(p) (4) 
v = A r p. (5) 

The stable finite element discretization for the PFC 
model in the plane introduced in [18] can be adapted 
to solve Eqs. ([3|-([5| on a surface triangulation using 
parametric finite elements. The key idea is to use the 
surface operators on the discrete surface which consists 
of triangles T. To do the integration on these triangles 
a parametrization Ft ' T — >• T is used, with T as the 
standard element in M 2 . These allow us to transform all 
integrations onto the standard element using the finite 
element basis defined also only in R 2 . The parametriza- 
tion Ft is given by the coordinates of the surface mesh 
elements and provides the only difference between solving 
equations on surfaces and on planar domains. For a sur- 
face we have to allow Ft : , whereas for a planar 
domain Ft : M 2 — >> R 2 . With this tiny modification any 
code to solve partial differential equations on Cartesian 
grids can be used to solve the same problem on a surface, 
providing a surface triangulation is given. The computa- 
tional cost is the same as solving the problem in a planar 
domain. Within an efficient implementation this does 
not even require to recompile a running two-dimensional 
(2D) code, but only a change in a parameter file, as e.g. 
done in AMDiS [19]. With this approach all available 
tools to solve partial differential equations on planar do- 
mains, such as adaptive refinement, multigrid algorithms 
or parallelization approaches are available also to solve 
equations on surfaces. 

The approach is used to evolve a randomly perturbed 
constant initial configuration p = po towards an energy 
minimum. With the possibilities to use adaptive time 
stepping and the efficiency of parallelization of the finite 




FIG. 1: (Color online) Local minimal energy configuration for 
6.064 particles on a sphere, (left) density profile, (right) color 
coded number of neighbors, 5-black, 6-red (gray) and 7- yellow 
(white). 



element method problem sizes of 1 x 10 6 particles can be 
addressed. The simulation for Fig. [l]with 6.064 particles 
required 1 day computing time on a single processor. 

In order to validate our approach we compute minimal 
energy configurations for various numbers of N. We sys- 
tematically compute the minimal energy configuration 
for all N e [12,2790]. In our configuration N = 2790 
corresponds to r = 100. The numerical results indicate, 
that the obtained minimal energy configurations are only 
sensitive to the defined lattice spacing and insensitive to 
a large extent to the other parameters. This might ex- 
plain why triangular tessellations on spherical surfaces 
occur in very distinct occasions, for which the interac- 
tions involved may differ a lot. In the following we use 
po = —0.3 and e = 0.4, which together with the radius 
r of the sphere determines N. For all N the type and 
number of disclinations, as well as the computed energy 
is in agreement with known analytical or other numerical 
results. For N < 112 the configurations and energies co- 
incide with the known equilibrium values. The maximal 
deviation from the minimal energy for larger N is less 
than 0.1 %. Table [I] shows computed configurations for 



TABLE I: Comparison with known results for small N. In 
order to compute the energy we identify the position of the 
maxima in the denisty field and compute the Coulomb energy 
according to this positions. 
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FIG. 2: (Color online) Annihilation of dislocation by local 
rearrangment of 5-7 defect. The gray lines indicate the for- 
mation of a new neighboure at an intermediate state. 



selected N. For N > 360 we obtain additional defects in 
the ground state, which are pairs of fivefold and sevenfold 
coordinated particles (dislocations) and chains of alterna- 
tion fivefold and sevenfold coordinated particles (grain- 
boundary scars). Since dislocations have vanishing total 
disclination charge there can be an arbitrary number of 
them in any spherical lattice configuration without vi- 
olating the topological constraint on the total disclina- 
tion charge discussed above. The grain-boundary scars 
found are unlike any grain boundaries found in flat space 
as they terminate freely inside the crystal at both ends. 
Due to the importance of dislocations and grain bound- 
aries in bulk materials in determining material proper- 
ties similar roles can be expected on surfaces. In [13] the 
motion of dislocations is observed experimentally. Dis- 
location motion can be separated into glide and climb, 
where glide is motion parallel to the dislocation's Burger 
vector and requires only local rearrangement of the lat- 
tice, and climb is motion perpendicular to the Burger's 
vector and requires the presence of vacancies and inter- 
stitials. In agreement with the observations in [13] and 
the computational results based on elasticity in [12] we 
also observe only glide motion which leads to significant 
shape changes of the scars. Fig. [2] shows a local rear- 
rangement of a dislocation. 

In [14] the question is asked which effect a raising of 
temperature has on the spherical crystal. For bulk poly- 
crystalline material there is indirect experimental evi- 
dence for the occurrence of grain boundary premelting, 
which could directly be visualized for colloidal crystals in 
[20] . As a liquid film at grain boundaries will alter macro- 
scopic properties and especially will lead to a drastically 
reduced resistance to shear stresses which can lead to ma- 
terial failure it is not only of theoretical interest if grain- 
boundary premelting is also present in spherical crystals. 
As grain-boundary scars belong to the equilibrium state 
it would be even more severe as it would be a general 
property of crystalline materials on curved surfaces. PFC 
models have already been used to study grain-boundary 
premelting in flat space [2TJ [22]. For high-angle grain 
boundaries an uniformly wetting is observed below the 
melting temperature in these studies. Raising the tem- 
perature (increasing e) but keeping it below the melting 
temperature according to the phase diagram indeed leads 
to melting of the grain-boundary scars, see Fig. [3] 

We use the premelting of the dislocations to improve 
the local minima our evolution settles in. Within an it- 




FIG. 3: (Color online) Time sequence showing premelting at 
grain boundary scars. The simulations indicate an initiation 
of the melting at the ends of the scars. 



erative procedure we evolve the system according to Eqs. 
([3| - ([5| until we reach a minimal state, increase the 
temperature and run the system until liquid layers have 
replaced all dislocations, and start the evolution again 
with the original temperature. The algorithm is termi- 
nated if the total energy does not further decrease. Typ- 
ically this is achieved after 2-3 iteration cycles. In Figure 
[4] we plot the excess dislocations in a scar as a function 
of the system size \ /f N = r/a. 
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FIG. 4: Excess dislocations as a function of system size. The 
obtained slope is 0.388 ±0.020, which is in excelent agreement 
with the experimental measured data in [4 , which give 0.404zb 
0.062, and the theoretical value of 0.41 from [10] 



As already pointed out the method is not restricted 
to spherical geometries. Indeed the algorithm works for 
arbitrary surfaces with the only requirement that an ap- 
propriate surface mesh is needed on which the computa- 
tion can be done. As an example we use toroidal crystals, 
which can be found e.g. in self-assembled monolayers of 
micelles and vesicles [24 or carbon nanotori [25 . We 
compute low energy configurations for various toroidal 
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lattice configurations and reach good comparison with 
results reported in [26] . Fig. [5] shows a typical configu- 
ration for aspect ratio R1/R2 = 2.78. 




FIG. 5: (Color online) Minimal energy configuration on torus 
with N = 239, there are 13 fivefold disclinations and 13 seven- 
fold disclinations, 5-black, 6-red (gray) and 7- yellow (white). 

In computations on more complicated surfaces, with 
convex and concave regions, we observe isolated fivefold 
and sevenfold disclinations which arrange according to 
the local curvature of the surface. Thereby fivefold discli- 
nations are preferably found in convex regions, whereas 
sevenfold disclinations are present in concave regions, 
which is in accordance with the theory discussed in [23] . 
The ability of the approach to work on arbitrary surfaces 
will also allow us to consider ordering on evolving sur- 
faces. As discussed above for low surface tensions a buck- 



ling transition of the disclinations can turn the sphere 
into a polygon inorder to reduce the energy. The ques- 
tion arises if such a transition can intervene with grain 
boundary scar formation. In [27 it is speculated that 
grain boundary scars could be formed on capsids at an 
intermediate stage of their evolution and that the release 
of the bending energy present in these scars into stretch- 
ing energy could allow for shape changes. To model such 
shape changes requires us to evolve the surface. Appro- 
priate continuum models which account for bending and 
surface tension are discussed in the mathematical review 
papers [28, 29 . Different concepts have been developed 
to solve differential equations on evolving surfaces, [30] in 
the context of parametric finite elements (as considered 
here) and [31] within a phase-field context. The coupling 
of the surface evolution with the evolution of the number 
density on it and with it the question if grain boundary 
scars can initiate a buckling transition, however remains 
open and requires further developments. 
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